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This paper describes a multidimensional hydro-dynamic code which can be used for 
the studies of relativistic astrophysical flows. The code solves the special relativistic 
hydrodynamic equations as a hyperbolic system of conservation laws based on High 
D ■ Resolution Shock Capturing (HRSC) Scheme. Two standard tests, one of which is 

n 

the relativistic blast wave tested in our previous paper and the other is the colli- 
sion of two ultrarelativistic blast waves tested in here, are presented to demonstrate 



^ • that the code captures correctly and gives solution in the discontinuities, accurately. 

m 

. The relativistic astrophysical jet is modeled for the ultrarelativistic flow case. 

\ The dynamics of jet flowing is then determined by the ambient parameters such 



O 



as densities, and velocities of the jets and the momentum impulse applied to the 



computational surface. We obtain solutions for the jet structure, propagation of jet 
$h . during the time evolution, and variation in the Mach number on the computational 

domain at a fixed time. 

X 

I. INTRODUCTION 

Many high-energy astrophysical problems (SNRs, 7-ray bursts -GRBs, Active Galactic 
Nuclei (AGN) hot spots, etc.) involve relativistic flows, and thus understanding relativistic 
flows is important for interpreting the astrophysical phenomena correctly. For instance, 
intrinsic beam velocities typically larger than 0.9c are required to explain the apparent 
superluminal motions observed in relativistic jets in microquasars in Galaxies^ as well 
as in extragalactic radio sources associated with AGN. The shocks created during these 
phenomena accelerate particles which emit the observed radiation. In particular, it is widely 
accepted that the recently discovered GRBs afterglow results from an emission by relativistic 
shocks, created by the interaction between an initial ejecta and the interstellar medium. 
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The recent observations of GRB afterglow have lead to numerous attempts to model these 
phenomena. 

General relativistic effects must be considered when strong gravitational fields are en- 
countered as, for example, in the case of coalescing neutron stars or near black holes. The 
significant gravitational wave signal produced by some of these phenomena can also only be 
understood in the framework of the theory of general relativity. Another field of research, 
where special relativistic flows are encountered, is heavy-ion collision experiments performed 
with large particle accelerators. The heavy-ions are accelerated up to ultra-relativistic veloc- 
ities to study various aspects of heavy- ion collision physics (e.g., multi-particle production, 
the occurrence of nuclear shock waves, collective flow phenomena, or dissipative processes) 
to explore the equation of state for hot dense nuclear matter, and to find evidence for the 
existence of the quark-gluon plasma. 

Multi-wavelength observation of extragalactic jets performed in radio, optical and X-ray 
bands have extended dramatically our knowledge about the complex phenomenology of these 
objects. By means of these observations, different kinds of objects that can house these jets 
have been characterized, at the moment we can classify distinctly these sources as radio 
galaxies and quasars. A deeper analysis of these and other related objects led us to the idea 
of AGN. The images obtained from Hubble Space Telescope revealed a new detail in gas 
flow and shock wave patterns involving astrophysical jets and colliding interstellar winds of 
particles. Astrophysical jets are defined as highly collimated outflows in the form of high 
velocity mass flows. The outflows are observed in young stellar objects, proto-planetary 
nebula, compact objects, AGN, and GRBs 

SB 

. The jets include high Mach number and 

interact with surrounding ambient gasQ|. 

Simulating the fluid flow and shock wave patterns and detailed temperature profiles 
by implementing theoretical models in a gas dynamics simulator will help in analyzing 
the processes at work in these astrophysical objects. In these investigations, we apply 
the HRSC scheme with appropriate initial and boundary conditions to simulate relativistic 
astrophysical jets with high Mach number from the compact objects. Pioneered in this 
field Ref. 0] is the pioneer in this field. They were able to show that a flow of supersonic 
plasma remains stable and develops features that could be identified with the features in the 
observations of radio galaxies. A circular and canonical deceleration area called Mach disk 
as the hot spot, which is a strongly collimated beam as the elongated structure of the jets, 



and a big zone of exhaust material as the lobes of radio galaxies Meanwhile, more physics 
has been included in the calculations such as the parameters varying in the perpendicular 
to the jet flow direction, the jet radius magnetic field 33 and the special relativistic 
effects llOlllli. Jets inside massive star, observed and called highly relativistic iets, have 
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been studied numerically in both Newtonian [12(1 [13J and relativistic [1J] simulations 
and it has been shown that the collapsed model is able to explain many of the observed 
characteristic of GRBs. The collapsar is formed when the iron core of a rotating massive 
star collapses to black hole and an accretion disk. The explosive deaths of massive stars 
produce relativistic jets. These studies are also declared that the interaction of jet with the 
matter at the stellar surface and the stellar wind could carry the information for activities. 
The cocoon of the jet would also have different properties and the shocks within cocoon and 
jet could lead to 7— ray and x— ray transients. 

In this paper, our goal is to simulate the ultrarelativistic jet problem using the already 
existing code which solves fully special relativistic hydrodynamics equations with HRSC 
scheme. 

II. FORMULATION 

A general technical description of the Hydrodynamic equation is given by DonmezQ as 
they will be used in our code development and the analytic description of different problems. 
The GRH equations are written in the standard covariant form, consist of the local conser- 
vation laws of the stress-energy tensor T^ v and the matter current density J M : 



V , t T^ = 0, V^ M = 0. (1) 

Greek indices run from to 3, Latin indices from 1 to 3, and units in which the speed of 
light c = 1 are used. Vm stands for the covariant derivative with respect to the 4-metric of 
the underlying spacetime, g^ v . 

Defining the characteristic waves of the general relativistic hydrodynamical equations 
is not trivial with imperfect fluid stress-energy tensor. The viscosity and heat conduction 
effects are neglected. This defines the perfect fluid stress-energy tensor. This stress-energy 
tensor is used to derive the hydrodynamical equations. Using perfect fluid stress-energy 
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tensor, we can solve some problems which are solved by the Newtonian hydrodynamics 
with viscosity, such as those involving angular momentum transport and shock waves on an 
accretion disk, etc. Entropy for perfect fluid is conserved along the fluid lines. The stress 
energy tensor for a perfect fluid is given as 

T» v = phu>*u v + Pg^. (2) 

A perfect fluid is a fluid that moves through spacetime with a 4-velocity which may vary 
from event to event. It exhibits a density of mass p and isotropic pressure P in the rest 
frame of each fluid element, h is the specific enthalpy, defined as 

h = l + e + -. (3) 
P 

Here e is the specific internal energy. The equation of state might have the functional form 
P = P(p, e). The perfect gas equation of state, 

P = (r - l)pe, (4) 

is such a functional form. 

The conservation laws in the form given in Eq.Q are not suitable for the use in advanced 
numerical schemes in 2D. In order to carry out numerical hydrodynamic evolutions and to 
use HRSC methods, the hydrodynamic equations after the 2+1 split must be written as a 
hyperbolic system of first order flux conservative equations. The Eq.fllJ is written in terms of 
coordinate derivatives, using the coordinates (x° = t, x l ) where and rest of the paper 2 = 1,2 
and j = 1,2. Eq.Q is projected onto the basis {n M , (^r) M }, where n M is a unit timelike 
vector normal to a given hypersurface. After a straightforward calculation conservation form 
of Special Relativistic Hydrodynamical (SRH) equation can be written, 

d t U + d x F x + d y F y = 0, (5) 

where dt = d/dt and d x = d/dx x . This basic step serves to identify the set of unknowns, the 
vector of conserved quantities U, and their corresponding fluxes F(U). With the equations in 
conservation form, almost every high resolution method devised to solve hyperbolic systems 
of conservation laws can be extended to SRH. 
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The evolved state vector U consists of the conservative variables (D, S x , S y , r) which are 
conserved variables for density, momentum in x and y direction and energy respectively; in 
terms of the primitive variables (p, v x ,v y ,e), this becomes, 



D = jyWn S x = ^phW 2 v x , S y = ^phW\- r = ^{phW 2 -P-Wp) (6) 

Here 7 is the determinant of the 3-metric 7^', which is unitary matrix, v x and v y are the 
fluid 3-velocity, the lapse functions = 1 and the shift vector j3 = for SRH and W is the 
Lorentz factor, 



W = au l 



(7) 



The flux vectors F % are given by 



/ 



F' 



v*D 



v i S j + P5) 



\ 



(8) 



/ 



(v l T + V l P 

The spatial components of the 4-velocity u l are related to the 3-velocity by the following 
formula: u l = Wv l . The source vector 5 = for SRH case. 



A. Spectral Decomposition and Characteristic Fields 

The use of HRSC schemes requires the spectral decomposition of the Jacobian matrix 
of the system, dF l /dU. The spectral decomposition of the Jacobian matrices of the SRH 
equations with a general equation of state was reported here . It is displayed the spectral 
decomposition, valid for a generic spatial metric, in the x-direction, (dF x /dll); permutation 
of the indices yields the y direction. 

We started the solution by considering an equation of state in which the pressure P is a 
function of p and e, P = P(p, e). The relativistic speed of sound in the fluid C s is given by 



r2 _ dP 
s ~ dE 



3 - = l + 7-h' (9) 



where x — 9P/dp\ e , k = dP/de\ p , S is the entropy per particle, and E = p + pe is the total 
rest energy density. 
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A complete set of the right eigenvectors [fi\ and corresponding eigenvalues Aj along the 
x-direction obeys 



dF x 
dU 



[r*] = \ { [rj] , i = l, ...,4. 



(10) 



The solution contains triply degenerate eigenvalues, 



\1 II \ X' 

A l — A 2 — A 3 



(11) 



The other eigenvalues are given as 



\x — 

± 1 - v 2 C? 



v x (l~C 2 s )± 



VCKl - ^ 2 )[(l - v*Cl) - v x v x (l - c!j] . 



A set of linearly independent right eigenvectors spanning this space is given as 



(12) 



rT 



lhW(K- P C* 



~ > v x , fj. , 1 



hW(K - pC 2 s ) 



(13) 



and 



(14) 



ft 



l,hW(v x 



A 



, TTr hW(l - v x v 2 
hWVy, 



v x X- 



(15) 



1 - v x \± 

Here, the superscript T denotes the transpose. 

Now, the eigenvalues and eigenvectors are computed in the y-direction using the infor- 
mation in the x-direction in Donmez for general relativistic hydrodynamic equations . 
Let M 1 (t/) be the matrix of the right eigenvectors, that is, the matrix having as columns 
the right eigenvectors with the standard ordering([J, Q|), M 1 = (r_, fx, fz, f+). To obtain 
the spectral decomposition in the other spatial directions q — 2(= y), it is enough to take 
into account the following symmetry relations: 



1) The eigenvalues are easily obtained by carrying out in equations and (|12|) 
a simple substitution of indices x by q. 

2) The matrix M 1 in the x-direction is permuted according to M q = V q i(lSA l ), 
where the two operators V q \ act on their arguments by the following sequential 
operations: 

• to permute the second and (q + l)th rows, 

• to interchange indices 1 = x with q. 

From Donat et. al. the left eigenvectors in the ^-direction are : 



h-W 



W 
K- 1 



Wv 1 



Wv y 



-W 



L 0,2 



he 



V v v 



[1 - v x v x ) 
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, . h 2 



hWV^ x + (±#)4 5) 



r M (i - K.AI) + (2/C - l)Vl{W 2 v x ¥ - T xx v 3 



T xy (l - KA X ± ) + (2/C - l)Vg(W 2 v y £ x - T 



(18) 



1 - JC)[-jv x + V x ± {W 2 ti x - T xx )} - icw 2 vu x 
The variables used for the left eigenvectors in the x-direction are defined as follows: 



K — Ct 



k = k/ p 



CI = v x - VI , VI 



v x — A± 



Al = 



1 — v x v x 
1 - v x \ x ± 



1 - A x ± = v x Vl , A x ± -A x = v x {C x ± - C% t 



{CI - C x ± ) + {AlVl - A X ± V X ) = 



a x = h 3 w(]c-i)(ci-c x :)c 



£ x = 1 — v x v x 

Detailed informations related with numerical solution of GRH equations using high res- 
olution shock capturing schemes in 3D are given in our paper Donmez^]. This reference 
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explains detail description of numerical solutions of SRH equations, Adaptive-Mesh Refine- 
ment, high resolution shock capturing method used, and solution of GRH equation using 
Schwarzschild coordinate as a source term. The code is constructed for general spacetime 
metric with lapse function and shift vector and is fully parallelized to make optimum use 
of supercomputers which is necessary to achieve the numerical modeling of real astrophys- 
ical problems, such as astrophysical jet problems, coalescing of the compact binaries and 
accretion disk around the compact objects. 

III. NUMERICAL RESULTS 

A. Boundary Condition 

Boundary conditions are set by filling the data in guard cells with appropriate values. In 
the numerical calculation, boundary filling plays an important role in the simulations. The 
computational grid is extended at both sides of the physical domain to compute the fluxes 
at interfaces. These extra cells are also called guard cells or ghost zones. In this paper, 
we have used outflow boundary condition at each ghost zone. This boundary condition has 
to be provided on each time step for all primitive and conservative variables in the special 
relativistic hydro code. The outflow boundary condition for the computational domain is 
as follows. For velocity: U£ = -\{\U?\ - C/f) and U^ +1 = + E7&). For the other 

variables: Uq = and U^ +1 = U^j, where n represents the time step. 

B. Relativistic Blast Wave 

Riemann problems with large initial pressure jumps produce blast waves with dense 
shells of material propagating at relativistic speeds. The Riemann shock tube is a useful 
test problem because it has an exact time- dependent solution, and tests the ability of the 
code to evolve both smooth and discontinuous flows. For appropriate initial conditions, 
both the speed of the leading shock front and the velocity of the shell material approach 
the speed of light producing very narrow structures. The accurate description of these thin, 
relativistic shells involving large density contrasts is a challenge for any numerical code. 
Some particular blast wave problems have become standard numerical tests. One of the 
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common test problem is called the special relativistic shock tube problem. Testing the code 
with this problem is given our paper DonmezQ). 

C. Collision of two relativistic blast waves in 2D 

The collision of two strong blast waves in one-dimension was used by Woodward et al. Q| 
to compare the performance of several numerical methods in classical hydrodynamics. In 
the relativistic case, Yang et al. Q| considered the same problem to test the high-order 
extensions of the relativistic beam scheme, whereas Marti et al. 2(| used it to evaluate the 
performance of their relativistic PPM code. In the last treatment, the original boundary 
conditions were changed (from reflecting to outflow) to avoid the reflection and subsequent 
interaction of rarefaction waves allowing for a comparison with an analytical solution. The 
initial data corresponding to this test, consisting of three constant states with large pressure 
jumps at the discontinuities separates the states. 

In this report, the initial data corresponds to collision of two strong blast waves in 2D, 
consisting of three constant states with large pressure jumps at the discontinuities separating 
the states (at r = 0.15 and r = 0.85, where r = sqrt(y 2 + z 2 )), as well as containing 
the properties of the blast waves created by the decay of the initial discontinuities which 
are located on the main diagonal of the square box and the data is extracted along the 
other diagonal axis. The initial data for this test is listed in Tabl^U The initial pressure 
discontinuities drive shocks into the middle part of the grid; behind them, rarefaction form 
and propagate toward the outher boundaries where they are falling out from computational 
domain. When the two shock waves collied at t = 0.41, a very dense shell is created. The 
collision gives rise to a narrow region of very high density as seen in Figffl bounded by two 
shocks moving at speeds 0.88 (shock at the left) and 0.711 (shock at the right). The velocity 
of shock wave during the time evolution is plotted in the right panel of FigJBl The presence 
of very narrow structures involving large density jumps requires very fine zoning and high 
resolution scheme to resolve the states properly. The high resolution shock capturing scheme 
we used here resolves the structure of the collision region satisfactorily well. FigJ2] shows 
that more resolutions are needed to resolve the preshock and postshock states reasonably 
well. The left panel at FigO] shows the behavior of pressure jump during the time evolution. 
The outflowing shock does not propagate along the grid and the contact discontinuity is not 
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created on the numerical results. 

D. Astrophysical Jet 

Investigation of gas flow in close binary systems is important in considering the accretion 
disk around the compact objects. The interaction of the disk with the compact objects, 
through processes that are still poorly understood, produces jets with high energy to mass 
ratio. Here we are primarily concerned with the propagation of these jets, not so much with 
how they are born. The propagating jets produce solar wind in the medium and this wind is 
highly supersonic. So obstacles in the flow, large variations in stream speeds, or fast ejecta 
produce shock waves, just as shock waves are produced by a supersonic jet. With a jet, 
sound waves develop ahead of the jet nose because the matter is traveling faster than the 
waves can escape. A shock wave forms, due to pressure buildup, to keep the air flowing 
around the jet. 

Jets are, by current theory, the best, perhaps, essential method of explaining how such 
large amounts of energy can be observed from objects at cosmological distances. However, 
there is much speculation as to the exact nature of these jets, the primary questions being: 
how do they arise and what is their structure. In this paper, we focus on the question, what 
is the structure of GRB jets? 

Here, we numerically model the 2D simulation of a relativistic jet propagating through 
an homogeneous atmosphere. The initial computational domain and values of beam for the 
flow velocity, Mach number, rest-mass density and adiabatic index are given in Fig|3] The 
fluid viscosity and thermal conductivity are neglected. The ambient medium is has a size 
of 10x120. The jet is injected at 4.8 < y < 5.2 in the direction of positive z— axis. We have 
constructed a slab jet, which is periodic in the x— direction. An ideal gas equation of state 
with adiabatic index 7 = 5/3 is assumed to describe both the jet matter and the ambient 
gas. In our calculation, 256 resolution in y and 4096 resolution in z are used for Figs J5l9l It 
is well known that the propagation of a supersonic jet is governed by the interaction of jet 
matter with ambient medium, which produces a bow shock in the ambient medium. The 
evolution of the jet was simulated up to the time when the head of the jet is about to reach 
to the boundary, which is far away from the injection point. 

The time evolution of relativistic jet for rest-mass density in the plane x = is given in 
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FigEl at different snapshots when the jet propagates through y — z plane. The presence of 
emitting matter moving at different velocities and orientations could lead to local variations 
of apparent relativistic motion within the jet. The distribution of apparent motions during 
the time evolutions is inhomogeneous and fringed as seen in FigE3 The structure of jet 
contains cocoon and vortexes with a bow shock. The cocoon contains shocked jet material 
deflected backward at the head of the jet. Since the jet is narrowly collimated, its beam is 
very thin. 

The FiglHl shows the variations in one-dimensional cut of rest mass density in the jet 
propagation direction, z, at final snapshot for two different jet velocities, v = 0.5 and 
v = 0.9. While the one with smaller velocity is called as a mildly relativistic jet, the 
other one is called the ultra-relativistic jet. The ultra-relativistic one displays rich internal 
structure with oblique shocks effectively decelerating the flow in the beam from Mach number 
equal to 140 at the injection point down to a value of about 10 around the head as seen in 
Fig0 The one-dimensional cut of the rest-mass density as given in FigJH] shows that the jet 
density propagates and oscillates on the working surface and it produces time varying waves 
and shock waves which interact with ambient medium. The pressure is larger at the head 
of jet and it causes the jet to propagate through the ambient medium at a larger speeds 
in both direction. The graphs seen in FiglHJshow the variation of the y velocity of cocoon, 
propagating through computational domain, at a fixed y = 5 along the z— axis for four 
different times. The head of jet has the largest y velocity which is represented with peaks 
as seen in FigJEl FigsEl El and U\ also depict that there is a Mach shock at the jet head, 
and material that has passed through the head is slowly forming a backflow along the jet, 
building up the a shear layer. 

We have also examined the propagation of two-dimensional relativistic jet in case of the 
different initial velocities and adiabatic values, T. The rest mass density for different T 
values are plotted in FigEl In agreement with previous studies on astrophysical problem, 
the relativistic jet is also more tightly packed with the smaller T. The lower T means cooler 
jet with larger Mach number of the flow. The cooler jet presents also a complex structure of 
internal shocks generated by pressure mismatches between the beam and the overpressured 
cocoon and by perturbations of the beam boundary by vortices and bulk motions within 
the cocoon. The cocoon is mainly formed by large vortices in cold jets with T = 5/3 while 
the strong beam collimation causes a large acceleration of the jet in cold jets with T = 4/3. 
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The beam gas is less efficiently redirected into the cocoon, and thinner cocoon with smaller 
vortices form with T = 4/3 as seen in FigEl 

IV. CONCLUSION 

We have described the main futures of stable 2D special relativistic hydrodynamical code 
and discussed test problems in Ref.Q and here involving strong shocks in two-dimension. 
The test problem solved here is called the collision of relativistic blast wave. The promising 
results from the test problem give us encouragement to simulate ultrarelativistic gas flows. 

We have carried out simulations of fully relativistic astrophysical jet in 2D. By detailed 
examination of a time series of ultrarelativistic flow, we have described the propagation of 
jet through the computational domain. The ultrarelativistic flow is thought to occur and to 
play a crucial role in the generation of GRBs. We have also computed the dynamic of jet 
in case of different adiabatic index used in equation of state. The jet with T = 4/3 form 
thinner cocoon with smaller vortices and the matter collimation causes a large acceleration 
of the jet. 
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TABLE I: Initial data for the two relativistic blast wave collision test problem in 2D. The decay of 
the initial discontinuities (at r = 0.15 and r = 0.85) produces two shock waves (velocities f s fc D cfc) 
moving in opposite directions. The gas is assumed to be ideal with an adiabatic index 7 = 5/3. 



Blast wave collision test problem 





Left 


Middle 


Right 


p 


1.0 


1.0 


1.0 


p 


100.0 


0.01 


10.0 


V x 


0.0 


0.0 


0.0 


Vy 


0.0 


0.0 


0.0 



16 




FIG. 1: The evolution of log of the density distribution for the colliding relativistic blast wave 
problem up to the interaction of the waves. The data is taken from the two-dimensional colliding 
wave extracting from the diagonal axis. The computation has been performed with relativistic 
hydrodynamical code on an equidistant grid of 1024 zones 
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FIG. 2: The evolution of density at t = 0.41 at different resolution. The preshock and postshock 
region are resolved with more grid zones. 




FIG. 3: The left pane: the evolution of pressure at two-different time. Right panel: the evolution 
of velocity of the blast wave. 
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FIG. 4: Initial representation of computational domain for astrophysical jet with the appropriate 
initial conditions. 
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FIG. 5: Contour maps of the logarithm of the rest-mass density for relativistic jet, moving with 
v = 0.9c, at different times. Snapshots of the proper rest-mass density distribution are plotted at 
y — z plane. 256 resolutions in y and 4096 resolutions in z are used. 
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FIG. 6: Vertical configuration values of rest-mass density along a cut to the z-axis, at r = 5.0, are 
shown for different jet velocity(mildly and ultrarelativistic cases). This figure illustrates the value 
of rest mass density during the evolution of jet flow at fixed time. Oscillatory jet rest mass density 
is achieved. 
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FIG. 7: Behavior of Mach number as a function of Z for different initial jet velocities at different 
times. It is taken at a fixed y=5.0 
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FIG. 9: 2D representation of jet depend of adiabatic index at a fixed time for different jet velocities 



